Partition Function Zeros and Finite Size Scaling 
of Helix-Coil Transitions in a Polypeptide 
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We report on multicanonical simulations of the helix-coil transition of a polypeptide. The nature 
of this transition was studied by calculating partition function zeros and the finite-size scaling of 
various quantities. New estimates for critical exponents are presented. 
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A common, ordered structure in proteins is the a-helix 
and it is conjectured that formation of a-helices is a 
key factor in the early stages of protein- folding It 
is long known that a-helices undergo a sharp transition 
towards a random coil state when the temperature is in- 
creased. The characteristics of this so-called helix-coil 
transition have been studied extensively [0, most re- 
cently in Refs. ||J4}. They are usually described in the 
framework of Zimm-Bragg-type theories || in which the 
homopolymers are approximated by a one-dimensional 
Ising model with the residues as "spins" taking values 
"helix" or "coil" , and solely local interactions. Hence, in 
such theories thermodynamic phase transitions are not 
possible. However, in preliminary work [0 it was shown 
that our all-atom model of poly-alanine exhibits a phase 
transition between the ordered helical state and the dis- 
ordered random-coil state. It was conjectured that this 
transition is due to long range interactions in our model 
and the fact that it is not one-dimensional: it is known 
that the ID Ising model with long-range interactions also 
exhibits a phase transition at finite T if the interactions 
decay like l/r CT with 1 < er < 2 ||. Our aim now is to 
investigate this transition in the frame work of a critical 
theory by means of the finite size scaling (FSS) analysis 
of partition function zeros. Analysis of partition function 
zeros is a well-known tool in the study of phase transi- 
tions, but was to our knowledge never used before to 
study biopolymers. 

For our project, the use of the multicanonical algo- 
rithm [Q was crucial. The various competing interactions 
within the polymer lead to an energy landscape charac- 
terized by a multitude of local minima. Hence, in the 
low-temperature region, canonical simulations will tend 
to get trapped in one of these minima and the simula- 
tion will not thermalize within the available CPU time. 
One standard way to overcome this problem is the ap- 
plication of the multicanonical algorithm |?J and other 
generalized-ensemble techniques (§] to the protein fold- 
ing problem H. For poly-alanine, both the failure of 



standard Monte Carlo techniques and the superior per- 
formance of the multicanonical algorithm are extensively 
documented in earlier work fl(][| . 

In the multicanonical algorithm (7) conformations with 
energy E are assigned a weight w mu (E) oc l/n(E). Here, 
n(E) is the density of states. A simulation with this 
weight will lead to a uniform distribution of energy: 



Pmu(E) oc n(E) w mu (E) = const 



(1) 



This is because the simulation generates a ID random 
walk in the energy, allowing itself to escape from any lo- 
cal minimum. Since a large range of energies are sampled, 
one can use the reweighting techniques |jll| to calculate 
thermodynamic quantities over a wide range of temper- 
atures by 



< A >T = 



dx A(x) w-^Eix)) e - 0E ^ 



dx w- 1 {E{x)) e- 0E{x) 



(2) 



where x stands for configurations. 

It follows from Eq. [I] that the multicanonical algorithm 
allows us to calculate estimates for the spectral density: 



n{E)=P mu (E)w~ I (E) 



(3) 



We can therefore construct the partition function from 
these estimates by 



Z(P)=J2n(E)u l 



(4) 



where u = e~' with (3 the inverse temperature, (3 = 
1/ksT. The complex solutions of the partition function 
determine the critical behavior of the model. They are 
the so-called Fisher zeros [ [H?|Jl3| ], and correspond to the 
complex extension of the temperature variable. 

Our investigation of the helix-coil transition for poly- 
alanine is based on a detailed, all-atom representation 



1 



of that homopolymer, and goes beyond the approxima- 
tions of the Zimm-Bragg model || . The interaction be- 
tween the atoms was described by a standard force field, 
ECEPP/2, 0] (as implemented in the KONF90 program 
p5[) and is given by: 



Etot — Ec + El j + Ehb + E tor , 
332q l q J 




Ui (1 ± cos(mxi)) ■ 



(5) 
(6) 

(7) 

(8) 
(9) 



Here, r*jj (in A) is the distance between the atoms i 
and j, and xi is the i-th. torsion angle. Note that with 
the electrostatic energy term Ec our model contains a 
long range interaction neglected in the Zimm-Bragg the- 
ory ||. Since one can avoid the complications of elec- 
trostatic and hydrogen-bond interactions of side chains 
with the solvent for alanine (a nonpolar amino acid), ex- 
plicit solvent molecules were neglected. Chains of up to 
N = 30 monomers were considered. We needed between 
40,000 sweeps (N = 10) and 500,000 sweeps (N = 30) 
for the weight factor calculations by the iterative proce- 
dure described in Refs. |t],[1(|. All thermodynamic quan- 
tities were estimated from one production run of N sw 
Monte Carlo sweeps starting from a random initial con- 
formation, i.e. without introducing any bias. We chose 
AT SW =400,000, 500,000, 1,000,000, and 3,000,000 sweeps 
for TV = 10, 15, 20, and 30, respectively. 

For our analysis of the partition function zeros we 
first divide the energy range into intervals of lengths 0.5 
kcal/mol. Equation |] becomes now a polynomial in the 
variable u and can be easily solved with MATHEMAT- 
ICA to obtain all complex zeros m° (j = 1,2,...). For the 
case of N = 10 we also repeated the calculation of the 
zeros for energy bin sizes 1.0 kcal/mol and 0.25 kcal/mol. 
The changes in the zeros were smaller than the statistical 
errors. The effect of the energy bin size on the zeros is 
also discussed in Ref. (H]]. Figure 1 shows the distribu- 
tion of the zeros for N — 30 and provides already strong 
evidence for a singularity on the real axis: in the case of 
the (analytic) Zimm-Bragg theory the zeros would be lo- 
cated solely on the negative real u-axis jl8| . We summa- 
rize in Table I the leading zeros for each of the four chain 
lengths, where we have used the mapping u = e _/3 / 2 due 
to our binning procedure. 

The FSS relation by Itzykson et al. G3] for the leading 
zero u\(N), 



AT 
iV 


Kc (it-, ) 


lax (iti ) 


Re (ft ) 


Im(ft ) 


10 


0.5620(60) 


0.0702(33) 


1.138(21) 


0.248(11) 


15 


0.6015(23) 


0.0472(21) 


1.0104(77) 


0.1566(67) 


20 


0.6105(29) 


0.03275(88) 


0.9842(94) 


0.1072(26) 


30 


0.6159(19) 


0.02200(78) 


0.9681(63) 


0.0714(25) 



TABLE I. First partition function zeros for poly-alanine 
chains of various chain lengths. 



shows that the distance from the closest zero it?, to the 
infinite-chain critical point u c — e~^°/ 2 on Re(it) axis, 
scales with a relevant linear length L, which we trans- 
lated as N 1 ^ in the above equation. Here, (3 C is the 
inverse critical temperature of the infinite long polymer 
chain and y is the correction to scaling exponent. We re- 
mark that, unlike in the Zimm-Bragg model, we have no 
theoretical indication to assume d as a particular integer 
geometrical dimension and report therefore estimates for 
the quantity dv. 

For sufficiently large N, the exponent dv can be ob- 
tained from the linear regression 



- \n\u\{N)-u c \ 



(11) 



This relation requires an accurate estimate for u c . There- 
fore, we prefer to calculate our estimates for dv from 
the corresponding relation with \u\ — u c \ replaced by its 
imaginary part ImuJ. Including chains of all lengths, 
N = 10 — 30, this approach leads to dv = 0.93(5), with 
a goodness of fit Q — 0.48. Figure 2 displays the corre- 
sponding fit. Omitting the smallest chain, i.e. restrict- 
ing the fit to the range N = 15 — 30, does not change 
the above result. We obtain now dv = 0.93(7), with 
Q = 0.22. This indicates that the dv determination is 
stable over the studied chains and therefore, the correc- 
tion exponent y can be disregarded in face of the present 
statistical error. 

Considering the real part of the leading zeros given in 
Table I, Re(/3°(iV)) = -^{[Reu^iV)] 2 + [imw°(/V)] 2 }, 
we can derive the critical temperature through the fol- 
lowing FSS fit pf, 



Re(/?°(A0) =l3 c + bN- 1 / d » . 



(12) 



We obtain (3 C = 0.906(12) (Q = 0.005) for the range N = 
10 - 30, and (3 C = 0.929(14) (Q = 0.63) for JV = 15 - 30. 
The last and more acceptable estimate, corresponds to 
T c (oo) = 541(8) K. 

A stronger version for the relation ([l0|) considers that 
the next zeros u°(iV), should also satisfy a scaling relation 



\u°JN) 



l/dv 



(13) 



+ AN- x ' du [l + 0{N y l d )) , y<0 (10) 



where j labels the zeros in order of increasing distance 
from u c . This relation is expected to be satisfied for large 



2 



N 


T 




r c 


Tmin 


b(T min ,N) 


10 


427(7) 


8.9(3) 


150(7) 


298 


-0.48(4) 


15 


492(5) 


12.3(4) 


119(5) 


429 


-0.59(10) 


20 


508(5) 


16.0(8) 


88(5) 


469 


-0.55(8) 


30 


518(7) 


22.8(1.2) 


58(4) 


500 


-0.20(4) 



TABLE II. Numerical results for poly-alanine chains of 
various lengths: critical temperature T c defined by the maxi- 
mum of specific heat C m ax , width Tc of peak in specific heat 
and temperature T m i n where the Binder cumulant b(T, N) 
has its minimum, b{T m i n , N). 



j and allows for an independent check of the estimate for 
our exponent dv. The scaling plot in Fig. 3 for the 
roots closest to the critical point u c demonstrates that 
the assumed scaling relation is indeed observed for our 
data as N increases and consistent with our estimate of 
the exponent dv. 

Our results for the critical temperature and critical 
exponent can be compared with independent estimates 
obtained from FSS of the specific heat: 



C N (T) = (3 2 (< E 2 {T) > -< E(T) > 2 )/N 



(14) 



Defining the critical temperature T C (N) as the position 
where the specific heat Cn(T) has its maximum, we 
can again calculate the critical temperature by means of 
Eq. [l^. With the values in Table II we obtain T c (oo) = 
544(12) K, which is consistent with the value obtained 
from the partition function zeros analysis, T c (oo) = 
541(8) K. Choosing T X (N) and T 2 (N) such that C(Ti) = 
1/2 C(T C ) = C(T<i), we have the following scaling relation 
for the width T C {N) of the specific heat pf, 



T C (N) = T 2 {N) - T X (N) oc N~ 1/du . 



(15) 



Using the above equation and the values given in Ta- 
ble II, we obtain dv = 0.98(11) (Q = 0.9) for chains 
of length N — 15 to N — 30, i.e. omitting the short- 
est chain. This value is in agreement with our estimate 
dv = 0.93(5), obtained from the partition function zero 
analysis. Including N = 10 leads to dv = 1.19(10), but 
with a less acceptable fit (Q = 0.1). The analysis of par- 
tition function zeros seems also to be more stable than 
one relying on Eq. |l^. No significant change in dv was ob- 
served when the data from ref. jij (which relied on much 
smaller number of Monte Carlo sweeps) were used in the 
partition function zeros analysis, while Eq. |l5| leads for 
this reduced statistics to an estimate for dv = 1.9. 

Through the scaling relation for the peak of the specific 
heat, we can evaluate yet another critical exponent, the 
specific heat exponent a, by: 



Cmax 
N 



X 



jya/dv 



(16) 



In particular, with the values for C^ ax as given in Ta- 
ble II, we obtain a — 0.86(10). The scaling plot for the 



specific heat is shown in Fig. 4: curves for all lengths of 
the poly-alanine chains nicely collapse on each other in- 
dicating the scaling of the specific heat and the reliability 
of our exponents. It worths to note that our estimates for 
dv and a, as obtained from the finite size scaling of the 
specific heat, obey within the errorbars the hypcrscaling 
relation dv = 2 — a. 

It is well-known that renormalization-group fixed point 
picture leads to a critical exponent dv = 1, a = 1 and 
7=1 for a first-order phase transition p9|-|2l|]. Our 
estimate dv — 0.93(5) for the correlation exponent devi- 
ates from unity and rather indicates that the 'helix-coil- 
transition" is a strong second order transition. However, 
the errorbars are such that a first order phase transi- 
tion cannot be excluded. Our values for the specific heat 
exponent a = 0.86(10) and the susceptibility exponent 
7 = 1.06(14) (data not shown) are consistent with a first 
order phase transition, but also not conclusive. A com- 
mon way to evaluate the order of a phase transition is by 
means of the Binder energy cumulant |22f , 



b(T, N) = 1 — 



< E A (T,N) > 



(17) 



3 < E 2 (T,N) > 2 ■ 

For a second order phase transition one would expect 
that the minimum of this quantity b(T m i n , oo) approaches 
2/3. Here T m i n defines the temperature where the cu- 
mulant reaches its minimum value and b(T m i n , oo) = 
limTv^oo b(T, N). With the present values of Table II 
we find the infinite volume extrapolation b(T m i n , oo) — 
0.23(13) (Q = 0.12), for the range N = 15 - 30, which is 
consistent with a first order phase transition. However, 
we cannot exclude the possibility of a second order phase 
transition because the energy cumulant scales with the 
maximum of specific heat Q, b(T, N) - jV«/*'-i j and 
the true asymptotic limit is reached only for rather large 
chains due to the value of a/dv. In fact, the straight 
line fit for the range N = 10 — 30 is less consistent with 
our data (Q ~ 0.001). Hence, we conclude that our re- 
sults seem to favor a (weak) first order phase transition, 
but are not precise enough to exclude the possibility of a 
second order phase transition. 

To summarize, we have used a common technique for 
investigation of phase transitions, analysis of the finite 
size scaling of partition function zeros, to evaluate the 
helix-coil transition in an all-atom model of poly-alanine. 
Although our results are due to the complexity of the 
simulated model not precise enough to determine the or- 
der of the phase transition, we have demonstrated that 
the transition can be described by a set of critical ex- 
ponents. Hence, we have shown for this example that 
structural transitions in biological molecules can be de- 
scribed within the frame work of a critical theory. 
Acknowledgements : 

Financial supports from FAPESP and a Research Excel- 
lence Fund of the State of Michigan are gratefully ac- 
knowledged. 



3 



Figure Captions: 



R.M. Ballew, J. Sabelko, and M. Gruebele, 
Proc. Nat. Acad. Sci. USA 93, 5759 (1996). 

D. Poland and H.A. Scheraga, Theory of Helix-Coil 
Transitions in Biopolymers (Academic Press, New York, 
1970). 

J. P. Kemp and Z.Y. Chen, Phys. Rev. Lett. 81, 3880 
(1998). 

U.H.E. Hansmann and Y. Okamoto, J. Chem. Phys. 110, 
1267 (1999); 111 1339(E) (1999). 

B.H. Zimm and J.K. Bragg, J. Chem. Phys. 31, 526 
(1959). 

F.J. Dyson, Commun. Math. Phys. 12, 212 (1969); 
J.F. Nagle and J.C. Bonner, J. Phys. C 3, 352 (1970); 

E. Bayong, H.T. Diep and V. Dotsenko, Phys. Rev. Lett. 
83, 14 (1999). 

B. A. Berg and T. Neuhaus, Phys. Lett. B 267, 249 
(1991). 

U.H.E. Hansmann and Y. Okamoto, in: Stauffer, D. (ed.) 
"Annual Reviews in Computational Physics VF (Singa- 
pore: World Scientific), p. 129. (1998). 
U.H.E. Hansmann and Y. Okamoto, J. Comp. Chem. 14, 
1333 (1993). 

Y. Okamoto and U.H.E. Hansmann, J. Phys. Chem. 99, 
11276 (1995). 

A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 
61, 2635 (1988); Phys. Rev. Lett. 63, 1658(E) (1989), 
and references given in the erratum. 
M.E. Fisher, in Lectures in Theoretical Physics, vol. 7c 
p.l (University of Colorado Press, Boulder, 1965). 

C. Itzykson, R.B. Pearson, and J.B. Zuber, Nucl. Phys. 
B 220 [FS8], 415 (1983). 

M.J. Sippl, G. Nemethy, and H.A. Scheraga, J. Phys. 
Chem. 88, 6231 (1984), and references therein. 
H. Kawai, Y. Okamoto, M. Fukugita, T. Nakazawa, and 
T. Kikuchi, Chem. Lett. 1991, 213 (1991); Y. Okamoto, 
M. Fukugita, T. Nakazawa, and H. Kawai, Protein Engi- 
neering 4, 639 (1991). 

U.H.E. Hansmann and Y. Okamoto, Physica A 212, 415 

(1994) . 

M. Karliner,S.R. Sharpe and Y.F. Chang, Nucl. Phys. 
B302, 204 (1988). 

V. Matveev and R. Shrock, Phys. Lett. A204, 353 

(1995) . 

M. Fukugita, H. Mino, M. Okawa and A. Ukawa, J. Stat. 
Phys. 59, 1397 (1990), and references given therein. 
M.E. Fisher and A.N. Berker, Phys. Rev. B 26, 2507 
(1982). 

N.A. Alves, B.A. Berg, and R. Villanova, Phys. Rev. B 
43, 5846 (1991). 

K. Binder, Phys. Rev. Lett. 47, 693 (1981). 

N.A. Alves, B.A. Berg, and S. Sanielevici, Phys. Lett. B 

241, 557 (1990). 



1 . Partition function zeros in the complex u plane for 
N — 30. For the Zimm-Bragg model the zeros 
would be located solely on the negative real u axis. 

2. Linear regression for -In (1m u\(N)) in the range 
N = 10-30. 

3. Scaling behavior of the first j complex zeros closest 
to u c = 0.6284, for chain lengths N = 10, 15, 20 
and 30. 

4. Scaling plot for the specific heat Cat(T) as a func- 
tion of temperature T, for poly-alanine molecules 
of chain lengths N = 10, 15, 20, and 30. 
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FIG. 1. Partition function zeros in the complex u plane for N = 30. For the Zimm-Bragg model the zeros would be 
located solely on the negative real u axis. 
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FIG. 3. Scaling behavior of the first j complex zeros closest to u c = 0.6284, for chain lengths N = 10, 15, 20 and 30. 
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FIG. 4. Scaling plot for the specific heat Cjv(T) as a function of temperature T, for poly-alanine molecules of chain 
lengths N = 10, 15, 20 and 30. 
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